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The symmetric two-layer Ising model (TLIM) is studied by the corner transfer matrix renormali- 
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be 1.773, which is consistent with the theoretical prediction 1.75 with 1.3% deviation. 
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a. Introduction The two-layer Ising model 
(TLIM), as a simple generalisation of the two- 
dimensional Ising model and a simple model for 
the magnetic ultra-thin films, has been studied for a 
long time §, |, |, |, |, |. Cobalt films grown on a 
Cu(100) crystal, for instance, have highly anisotropic 
magnetisation Q and could therefore be viewed as 
layered Ising-spin systems. It has been found that 
capping PtCo in TbFeCo to form a two-layer structure 
has applicable features, for instance, raising the Curie 
temperature and reducing the switching fields for 
over- writable magneto-optical disks ||. Apart from 
various possible applications to real physical materials, 
the model is theoretically interesting for its rich phase 
structure. The model has several interesting equivalents, 
such as a two-species gauge invariant Ising model [5J, a 
spin- 1 Ising model lief , a model of the dilute lattice gas 
and a quantum spin-i ladder at zero temperature. The 
TLIM is important for the investigation of crossover 
from the 2-dimensional Ising model to the 3-dimensional 
one. In particular, it has been argued that the critical 
point of the latter could be found from the spectrum of 



the TLIM |11 



In recent years, some approximation methods have 
been applied to this model |o|, || || ||, £jj[ 0, || . 
A critical line has been found in all these studies. As 
expected, the Curie temperature is very sensitive to the 
inter-layer interaction. Many discussions have been di- 
rected to the shift exponent at the decoupling point. Abe 
H and Suzuki j| have predicted 7 = 7/4 for the isotropic 
model many years ago. Only in recent years the compu- 
tational results have converged to the prediction with still 
about 2.3% deviation M, ifl. 



'Supported by the National Natural Science Foundation of China 
under the project 19772074 and by the Deutsche Forschungsgc- 
mcinschaft under th e project Schu 95/9-3. 
t Electronic address: LSchuclkc@t-onlinc.dc 



Apart from the shift exponent, it is also interesting to 
study the critical behaviour along the critical line. The 
model has the same critical exponents at the two ends of 
the critical line corresponding to the solvable decoupling 
limit and the infinite interlayer coupling limit. But it is 
clear that the decoupled system has a higher symmetry 
than the coupling layers, hence one cannot assume a pri- 
ori that the TLIM belongs to the one-layer Ising univer- 
sality class. Unusual finite-size effects have been observed 
in the dynamic process jl9). It has been proposed that 
the critical exponents (or some of them) would vary con- 
tinuously along the critical line ||, |l^]. However, there 
are also arguments in favour of unchanged exponents. 
Based on their correlation length equality method, An- 
gelini et al. argued that the exponent v is a constant [fl3| . 
The accurate computation of critical exponents along the 
critical line remains a difficult unsolved problem. 

It is our purpose to provide a reliable prediction for 
the critical exponents based on the corner transfer matrix 
renormalisation group (CTMRG) method |(| ||||. 

The CTMRG method follows White's idea of density 
matrix renormalisation group |23|, p4| . It reduces the 
phase-space dimensions by omitting the eigenstates cor- 
responding to small eigenvalues of the density matrix. 
But instead of constructing the density matrix through a 
series of column transfer matrices, Nishino and Okunishi 
made use of the corner transfer matrices. The CTMRG 
method has shown many merits for the two-dimensional 
classical models which can be mapped into interaction- 
around-face (IRF) models. Since the system can be eas- 
ily extended to very large lattices, the finite size effect 
can be ignored. The error mainly comes from the finite 
dimensions of the renormalised phase space, which only 
becomes severe for huge lattices at or very close to the 
critical point. Numerically it has been observed that the 
extrapolation from the finite size scaling can provide very 
accurate information about the criticality. 

We will focus on the TLIM composed of two identical 
layers. In this case, the model has a symmetry of inter- 
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changing the spins of the two layers. This symmetry is 
vital to the critical behaviour of the model. A breakdown 
of universality due to a small symmetry breaking flTj has 
been observed. In order to keep this symmetry, special 
care should be taken in the programming. Technically, 
the symmetry reduces the dimensions of the transfer ma- 
trix remarkably. But the extra degeneracy will slow down 
the convergence of the iteration which makes the prob- 
lem much more difficult compared to the one-layer Ising 
model. 

The critical behaviour can be easily observed in various 
quantities, such as the rapid increase of the magnetisa- 
tion, the lambda peak of the capacity, and the anoma- 
lous slowing down of convergence of the renormalisation 
group iterations. However, to calculate the critical ex- 
ponents, one needs to locate the critical points to high 
precision. This is done by a careful extrapolation of the 



renormalisation dimension m. 

In Section 2, the TLIM is defined and transformed into 
an IRF model which is suitable for the CTMRG study. 
In the same section, the CTMRG method particular for 
the considered model will be introduced. The finite-size 
scaling and the finite renormalisation dimension effect is 
studied in Section 3. At the same time the critical points 
are located, critical exponents of the scaling law are com- 
puted. Results for the critical exponents r/ and v can be 
found also in Section 3. Discussions and conclusions are 
given in the last section. 



b. The CTMRG method for the TLIM The 

Hamiltonian of the model is defined by 
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where Si = ±1 and er.; = ±1 are Ising spins on two iden- 
tical square lattices. For concreteness, let us say, Si is 
at the site i of layer 1, <ii is at the corresponding site of 
layer 2 which is also labelled by i. Then < i, j > should 
be understood as a nearest neighbour pair of sites either 
on layer 1 or layer 2, depending on which of the spin vari- 
able, s or (j, is concerned. In this paper, we only discuss 



the symmetric ferromagnetic case with nearest neighbour 
coupling Ji = J2 = J > and with interlayer coupling 
A = pJ > 0. 

To obtain the equivalent IRF model, two species of 
Ising spins, itjj = ±1 and — ±1, are introduced for 
each link < i, j >. The Hamiltonian (nf) is generalised to 
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-J' 
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E i <> 
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Vij ((Ji +(Jj) + p'y j s i a l 



(2) 



where J' and p' are defined as 

J' = ilii(e !J + \/e«-l), (3) 

It is not difficult to show that, by summing over all u— 
and v— spins, the TLIM is recovered apart from a overall 
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irrelevant constant factor to the partition function. In- 
versely, if the s-spins and cr-spins are summed, one will 
obtain the IRF model with four u-spin interaction and 
four v-spin interaction around each vertex. The vertex 
weight of the IRF model is given by 



W a 



bed 



^ exp { J' [s(u a + u b + u c + Ud) + (r(v a + v b + v c + v d ) + p'sa] } 



(5) 



8,<T=±1 



where a, 6, c, and d denote the four (u,v)-spin pairs in 
the links joining to the vertex. Special weights in the 



boundary, which have links less than four, can be defined 
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in a similar way as (|1) but the values of s and a are 
subject to the boundary condition. For instance, for the 
fixed boundary condition, the weight along the boundary, 
W abc , is defined by 

W abc = exp [J' (u a + u b + u c + v a + v b + v c + p')} . (6) 

The vertex at the corner has a weight 

W ab = exp [J' (u a + u b + v a + v b + p')] . (7) 

They are schematically given in Figure [l]a,b,c, respec- 
tively. 
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FIG. 1: (a) A full vertex weight W a bcd', (b) an edge vertex 
weigh W a bc', (c) a corner vertex weigh W ab - 

The partition function is a trace of products of all ver- 
tex weights. The trace here means to sum over all spins. 
In the following, whenever two vertices are connected by 
a link, a summation over the (u,v)-spin pair in that link 
is implied. The vertex weights of a quarter lattice, af- 
ter summation over all inner spins and boundary spins 
as it has been implied, is called a corner transfer ma- 
trix(CTM), denoted by C a p(N). Here a and j3 denote 
the spin configurations of two open sides of the corner re- 
spectively, and N is the size of the corner. A schematic 
representation for a CTM is given in Figure Ha. 
I 
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FIG. 2: (a) A corner transfer matrix (CTM), C a/3 (2); a 
half-row transfer matrix (HRTM), P a p, a (2). 



Another matrix we need is the so-called half-row trans- 
fer matrix (HRTM) denoted by P a p <a {N). It is defined 
by the recursion relation 



W abcd P a>l3> . c {N -l) 



(8) 



Ac,a(l) = W abc , (9) 
where a' and (3' are spin configurations of two sides of 
the HRTM with size of {N — 1) ; c is the (u,v)-spin pair of 
the top of that HRTM; a and (3 are spin configurations 
of the enlarged HRTM with a representing (a',b) and 
(3 representing (0',d) respectively. For each value of a, 
Pa/3,a(N) can also be considered as a matrix, denoted by 
P a (N). Figure |b is an example of the HRTM. 

In the process of the CTMRG iterations, the size of 
CTM is enlarged by one lattice spacing each time. This 
is done in this way: (1) glue a HRTM to each open side 
of the CTM, (2) complement the bigger corner by adding 
a vertex weight. The recursion relation is 



a" /3" cd 



P a >0»AN - 1) .C a „p„{N - l)P a »p>, d ( N - !) 



(10) 
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It is schematically shown in Figure || 

By use of CTM, the partition function of an even size 
lattice is simply the trace of the CTM to power four 



Zin = Tr(C{Nf). 



(11) 



Therefore, the CTM and the density matrix have com- 
mon eigenvectors. Denote the eigenvalues of C by 
{ujk\k — 0,1,2,...}, assuming that ojq > lo\ > LO2 > 
the eigenvalues of the density matrix are {^^Ifc = 
0,1,2,...}. According the spirit of DMRG, only a certain 
number of the biggest eigenstates will be kept. For in- 
stance, the first to eigenstates are kept, the eigenstates 



with smaller eigenvalues {o;fe|/c > m} are ignored, m is 
the so-called renormalisation dimension. 

The CTMRG iterations contains the following steps: 
(1) construct C(N) of a small CTM which is small enough 
to be exactly diagonalised; (2) construct P a (N); (3) di- 
agonalise C(N); (4) if the dimension of C(N) is smaller 
than m, all eigenstates are used as basis of the renor- 
malised phase space, otherwise only to eigenstates with 
the biggest eigenvalues are used as basis; (5) by use of 
the eigenvectors which have been chosen for the renor- 
malised space, form a projection operator U(N) that 
projects old spin states into the renormalised space; (6) 
use U(N) to project P a (N) to the renormalised HRTM 
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FIG. 3: Enlarging C Q/3 (2) into C a g(3). 



P£(N) ; (7) form a diagonal matrix C r (N) by the first m 
largest eigenvalues, which is the renormalised CTM; (8) 
use P£(N), C r (N), and the vertex weight W to form a 
bigger CTM, C(N + 1), and a bigger HRTM, P(N + 1), 
and repeat the above procedure from step (3) with N 
replaced by N + 1 . 

Since the renormalised phase space has dimension < 
m, the CTMRG in principle can infinitely go on. There- 
fore the lattice can be easily enlarged to very big sizes. 



Particularly for the TLIM, it is efficient to choose the 
basis to be symmetric and antisymmetric in the exchange 
of two layers. In this way, the CTM is automatically 
block diagonal. The symmetric block and the antisym- 
metric block can be diagonalised separately. This is not 
just allowing larger m and high precision, but is also es- 
sential for keeping the presumed symmetry between two 
layers. 

For small interlayer coupling, there are many approx- 
imate degenerate eigenstates due to the layer symmetry. 
If one does not separately consider the symmetric and 
antisymmetric phase space, it is easy to destroy the layer 
symmetry in the renormalisation procedure. 

Following the method of Nishino and Okunishi ptj , 
the magnetisation M can be calculated on an odd size 
lattice which is constructed by inserting a HRTM matrix 
between two adjacent CTM matrices and adding a proper 
vertex at the lattice centre. The probability that the 
central site has a spin state (s, a) is given by 



M s 



Eg bcd KlcdTr (P^C r P b r C r P^C r P^Cn 

Eabcd WabcdTr {PrCrprCrprCrpr C r) 



(12) 



where X*% cd is defined by 



x lbcd = ex P {J' i s ( u a + u b + u c + u d ) + a(v a + v b + v c + v d ) + p'sa]} 



(13) 
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The magnetisation of layer 1, M =< s >, is 

M = M++ + M+- - M - M-"* 



(14) 



Due to the layer symmetry, it is also one half of the to- 
tal magnetisation. The local energy density E can be 
calculated in a similar way. 



E 



where 



Eabcd YabcdTr (prC r P£C r PrC r PSCn 

Eabcd WabcdTr (P^C r PI C r Pi C r P r d C r ) 
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^ Xabcd 
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c. Extrapolation and finite-size scaling At the 

critical point J c , it is expected that M obeys the finite- 
size scaling form M ~ L~( d ~ 2+V ^ 2 , where L is the lattice 
size, so 



for the critical point J c with the vertical value of the 
line as — (d — 2 + rf)/2. For J near J c , there should be 
deviations to the power-law behaviour. A simple physics 
consideration reveals that the curves should be in the 



J=Jc 



InL 



FIG. 4: Schematic plot of the ideal curves din M /d\n L 
versus In L. 



\txM{L) — const — 



d-2 



n 



(17) 



If we calculate <91nM/<91nL and plot it with InL as 
the horizontal axis, the curve should be a horizontal line 



form as shown in Figure |J. So by searching for the "best" 
horizontal line, one can determine the critical point J c as 
well as the critical exponent rj. 

Things are not so easy, however. Figure |5| shows the 
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4 5 6 7 8 



In L 

FIG. 5: d\nM/dhxL plotted versus InL for p = 0.1, J = 
0.397726. The solid curves are for m = 60, 80, 120, 240. The 
dashed- line curve is the result of an extrapolation to m — *• oo. 

curves of dhxM / <91nL against InL for p = 0.1, J = 
0.397726 with different renormalisation dimensions m. 
This J value is quite near to J c . For small L, there exists 
a region where the lattice size is too small and finite-size 
scaling has not come into play. We will refer to this re- 
gion as the small-size region. As L increases, the number 
of states of the system increases while the renormalisa- 
tion dimension m remains unchanged, and we will suffer 
from the "finite-m effect" in due time. To get a more 
complete knowledge of the finite- m effect, let us look at 
Figure where curves for different J with m — 60 are 
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FIG. 6: dlnM jd InL plotted versus InL for p = 0.1, m = 60, 
for three different values of J. 

shown, m — 60 is relatively small in our calculation and 
it has a more drastic finite-m effect. When L is small, all 
the three curves shown in Figure ^ increase, which is an 
indication of the small-size effect. After a certain value of 
L, one of the curves begins to bend down and the other 
two become much flatter. Here the critical behaviour 
dominates. However, the finite-m effect soon shows itself 
by raising the curves, which is also demonstrated clearly 
in Figure [5]. As a result, only approximately the part in- 
side the dotted-line box resembles Figure ^. It is hard to 
tell exactly where the finite-m effect begin to dominate 
without comparison with curves of higher m. 



m 


Ji m) (d 


-2 + ij)/2 


60 


0.3977212 


0.12629 


80 


0.3977230 


0.12620 


120 


0.3977252 


0.12548 


240 


0.3977260 


0.12529 


oo 


0.3977264 


0.12520 



TABLE I: The critical point and critical exponent (d— 2+r/)/2 
at p — 0.1 estimated through first looking for the most flat 
curve of d In M/d In L ~ In L for certain m then extrapolating 
to infinite m. 



Nevertheless, some useful information can be drawn 
from the analysis. If a curve bends down for a period of 
L, one can judge that it corresponds to a J that is smaller 
than J c . In other words, if a curve has a local maximum, 
it corresponds to a J with J < J c . Thus the point now 
is to find Jc 60 ^ , above which the corresponding curve has 
no local maximum and below which it has. ji 60 ^ can 
be understood as a lower limit for J c . In general, one 
can obtain Jc™^ as an estimate of J c . The former will 
approach the latter from below as m — * oo. Meanwhile 
the vertical value of the horizontal part of the curve with 
J = Jc" 1 "* gives the estimate for (d — 2 + r))/2. To make 
the above discussion simpler, we have limited it only to 
the case with similar behaviour as that of p = 0.1. Later 
we will see in Figure that the small-size effect may 
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FIG. 7: 91n Af/dlnL plotted versus InL for J — J c with 
p — 1.0, 0.4, 0.1 from above, respectively. 

show itself differently as p changes. But similar analyses 
can still be developed for other cases. 

In most of our calculations, we chose m to be 60, 80, 
120, and 240. For each m, we calculated the magneti- 
sation with a set of J values and then determined the 
critical point using the method described above. The re- 
sults for p = 0.1 are shown in Table |. The last row shows 
the extrapolated values for m = oo. 

An alternate and finer way for dealing with the finitc- 
m problem is to extrapolate the data to m — ^ oo first 
and then analyse the extrapolated data. The extrapo- 
lation should lead to reasonable and smooth curves for 
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m = oo, and we assume the data converge monotonically 
as m — » oo. This defines our criterion for choosing the 
extrapolation formula. We still take as an example the 
case with J = 0.397726, p = 0.1, whose curves are plot- 
ted in Figure ||. For each value of lattice size L, we fit 
the points (l/m,M^ (£)), where m = 60,80,120,240, 
with 

y = a^x° + a^x A + a^x 3 + clq . (18) 

ao then is considered to be the value for m = oo (see 
Figure ^). Eq. (|i~8|) does not include terms of x 1 and x 2 




1/DimRen 

FIG. 8: Extrapolation from m = 60,80,120,240 for J = 
0.397726, p = 0.1 and L = 2001. 

because including them violates the criteria above. The 
extrapolated curve is also shown in Figure ^| as a dashed 
line. 

Now that we obtained the curves for m — oo with 
different J, we can search for the J value whose 
dlnM/dlnL ~ lnL curve is closest to a horizontal line. 
This J value is then the estimate of J c . In our calcu- 
lation, we use interpolation technique to help locate the 
"best" J, as shown in Figure ||. The solid lines are those 
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FIG. 9: dlnM/dlnL plotted versus InL for p = 0.1. The 
solid lines are the extropolated curves for J = 0.397726 and 
J = 0.397727. The dashed line is obtained from an interpo- 
lation between the two values of J. 

extrapolated to infinite m. We can obtain the magnetisa- 
tion for any value of J that is between J — 0.397726, for 
which the curve bends down, and J — 0.397727, where 



p 


(d 


-2 + r,)/2 


d — - 


0.1 


0.3977262 


0.12507 


1.0014 


0.4 


0.3541412 


0.12503 


1.0001 


1.0 


0.3117577 


0.12502 


0.9996 



TABLE II: Critical point J c , critical exponents for p — 
0.1,0.4, and 1.0, estimated through first extrapolating the 
curve of 91nM/<91nL ~ InL for infinite m then looking for 
the most flat curve. 



the line bends up, by linear interpolation. We found J c 
to be 0.3977262 with (d - 2 + r?)/2 = 0.12507 . 

The small-size effect for different ratios is presented in 
Figure [?]. For small p, this effect is big. In the case of 
p = 0.1, for instance, the small-size effect is visible with 
L up to 800. However, it seems that we overcome the 
small-size effect with L — 2000. 

Having J c in hand, we can determine another critical 
exponent v through the finite-size scaling form of the 
energy density 

E{L) - ~ L-( d -^ . (19) 

The critical point J c , critical exponents (d — 2 + rj)/2 
, and d — ^ for p — 0.1,0.4, and 1.0 are collected in 
Table |n[ 

We also estimated the shift exponent by taking p = 
0.005 and 0.01. The critical points are found to be 
J c = 0.432507 and 0.428593 respectively, which lead to a 
value 1.773 for the shift exponent that consists with the 
theoretical prediction 1.75 with 1.3% deviation. 



d. Discussion and conclusions Our numerical 
results suggest that the TLIM model is in the Ising uni- 
versality class. That is, for a finite inter layer coupling, 
the critical behaviour of the TLIM is the same as that of 
the one-layer system. Both < s > and < a > are order 
parameters in the transition. Our observations strongly 
exclude a phase transition for the electric order, except 
at the decoupling point p = 0. 

The critical line of the TLIM is in fact not a fixed 
line since the interlayer coupling is a relevant interaction 
which should flow to infinity in renormalisation group 
transformations. Then the TLIM resembles the one-layer 
Ising model in the long-range critical behaviour. To con- 
firm this picture, we also carry out the real-space renor- 
malizational group transformation to the model. The 
lattice is divided into 3x3 blocks. All nearest neighbour 
interactions are kept. The RG trajactories are plotted in 
Figure [l^. Two attractive fix points both having infinite 
p are observed. They are the ordered and disordered fix 
point, respectively. The critical line seperates two attrac- 
tive fix points and obviously it is not a fix line. 

However, as the interlayer coupling is small, it needs 
a large scale transformation to reach the infinite fixed 
points, so the model has extraordinary large finite-size 
effects as we have observed in the numerical results. 
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FIG. 10: The RG trajectories mapped to the TLIM param- 
eter space are plotted. The vertical axis is exp( — 1/A) with 
A = pj. Corresponding to the curves in the order of the leg- 
end, the starting couplings are J = 0.477, 0.438, 0.438, 0.439, 
0.42, and 0.35, respectively. 



istence of the severe competition between the unstable 
decoupling fixed point and the stable infinite interlayer 
coupling fixed points which, however, are far away. This 
supported in some sense by the unusual sensitivity of the 
shift exponent on the couplings Jl6|, [llj and the strange 
layer symmetry breaking in the effective state at the de- 
coupling limit Jl4| . Our estimate for the shift exponent 
is consistent with the theoretical prediction. 

The unusual small-size effects will have important con- 
sequences on the local properties of the TLIM. It is 
manifested in the short-time critical behaviour that we 
recently observed in Monte Carlo simulations which will 
be reported elsewhere [tjjl. 
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